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ABSTRACT 

Although different approaches to model a polarimctcr's accuracy have been described before, a complete error 
budgeting tool for polarimetric systems has not been yet developed. Based on the framework introduced by 
Keller & Snik, in 2009, we have developed the M&m's code as a first attempt to obtain a generic tool to 
model the performance and accuracy of a given polarimeter, including all the potential error contributions and 
their dependencies on physical parameters. The main goal of the code is to provide insight on the combined 
influence of many polarization errors on the accuracy of any polarimetric instrument. In this work we present the 
mathematics and physics based on which the code is developed as well as its general structure and operational 
scheme. Discussion of the advantages of the M&m's approach to error budgeting and polarimetric performance 
simulation is carried out and a brief outlook of further development of the code is also given. 
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1. INTRODUCTION 

Polarimetry is a very valuable remote-sensing technique that often yields information that is unobtainable through 
other techniques and is used in many different fields such as astronomy, Earth observation, biomedical diagnosis 
or land-mine and target detection, among others. 

Nevertheless, the construction of such instruments has often relied on the knowledge obtained through experi- 
ence than on a formal systems engineering approach, which is common practice when designing optical (imaging 
or spectroscopic) systems. Also, quite often, the polarimetry is implemented as an add on to an existing system 
which makes it sub-optimal by definition. It is only now, due to the increasing size and complexity of current in- 
strumentation projects, that this is starting to be demanded. For this we need to be able to predict the behavior 
and accuracy of the polarimetric system in order to optimize the design process. Only with the implementation 
of polarimetric error budgeting can the polarimetric performance be fully traded off against optical performance 
and other merit functions (like cost). 

In the systems engineering approach, the design process starts by setting the scientific requirements of the 
instrument which then must be translated into accuracy and sensitivity requirements. Polarimetric sensitivity 
is defined as the smallest signal that the instrument can detect above the noise and polarimetric accuracy is the 
uncertainty with which the instrument is able to measure a polarization signal after it has been detected with 
sufficient S/NP^l The former is limited by (photon) noise and spurious polarization effects that can be created 
by e.g. variable atmospheric properties or source variability. These effects are to a large extent decoupled from 
the polarimetric accuracy, and formalisms and codes exist to model themP The accuracy is limited by "real" 
polarization errors, which can be described by Mueller matrices. The goal of the design process then is to assure 
that the instrument will indeed meet the accuracy requirements when operating. 

With this aim, one needs to simulate the performance of the different preliminary designs and make an 
estimation of the error budget throughout the system to detect which elements contribute, and in which way, to 
the final response of the instrument. This allows to deal with the limitations of the system in an early design 
stage. In optical systems dealing only with intensity it is common practice to make use of error budgeting and 
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performance simulation tools that can carry out such an optimization process. Often this optical error budgeting 
involves adding wavefront errors (which are small, independent and scalar) in an RSS (root sum square) fashion, 
and force the total error to be smaller than a certain required value. The individual errors can then be distributed 
top-down, or added up bottom-up by adopting measured or modeled values, or a combination of both. So far, to 
our knowledge, a complete polarimetric error budgeting tool, working with a library of all possible error sources 
and their possible interactions has not yet been developed for polarimetric systems and this is mainly due to 
the complexity of the error propagation in these systems. Tyo (2002^1 an( j Boger et al. (2003^1 introduced 
polarimetric error budgeting with similar mathematical formalisms as presented in this paper, but their scopes 
focused on specific polarimetric elements. 

Errors in polarimetry have to be expressed as vectors and their values are often larger than the measured 
signal itself, due to the fact that some elements, essential to perform polarimetry, affect the polarization state of 
the incoming light in a major way. Often the degree of polarization of the signals to be measured is very low, 
and possibly much lower than the instrumental polarization. The main implication of this is that the common 
algorithms applied to optical systems for error propagation as RSS cannot be applied as such to polarimetric 
systems. In 2009 Keller & Snik analyzed this problem and developed a mathematical framework to transform 
the errors into additive ones to make them suitable for error budgeting and estimate the contribution of each 
physical parameter to the overall matrix.^ 

Based on this framework we have developed the M&m's code that computes the error propagation through the 
polarimetric system of all the potential error contributions, e.g. misalignment and varying material properties, 
and their dependencies on (global) physical parameters, e.g. wavelength and temperature. The ultimate goal 
of the code is to provide insight of the contribution of error sources to the final polarimetric performance and 
estimate the polarimetric accuracy of a given design. The code only pertains to predicting the polarimetric 
accuracy of a certain instrument as it relates the incoming Stokes vector to the measurement result. 

In Section [2] a description of the mathematical framework and the computational procedure is given. In 
Section [3] the physics for the modeling of error sources that M&m's uses in its library are explained. In Section 
[4] a description of the structure and operational modes of the code is given. Section [5] discusses the advantages 
and disadvantages of this approach for error budgeting of polarimetric systems and the next steps to be taken 
in the development of the simulator. 

2. ERROR PROPAGATION IN POLARIMETRIC SYSTEMS: THE MATH 
2.1 Mathematical approach 

In the Stokes formalism, every element in an optical system can be described as a 4 x 4 matrix called Mueller 
matrix (M e i ement ). The Stokes vector going out of each element S out = (I ut, Qout, U out , V out ) T is then obtained 
by multiplying this matrix by the incoming Stokes vector Si n = (Ji n , Qi n , Ui n , Vi n ) T , viz. : 

Sout = M-elementSin • (1) 

Each element of this matrix represents a relation between the components of S; n and S out : 



M 



element 



^ ^ ^out Qin ^ 3-out U in V lout Vi n ^ ^out 

^in ^ Qout Qin ^ Qout Ui n Y Qout Vin ^ Qout 
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(2) 



If we have a set of n elements forming an optical system we then have to multiply the respective element 
matrices from the last to the first optical element on the light path to get the total matrix for the system M to t ■ 

M^t = M n M n _ 1 ...M 2 M 1 . (3) 



Which obviously relates the incoming and outgoing Stokes vectors as follows: 



S ut — MtotSin ■ 



(4) 



Each element's performance, and hence its corresponding Mueller matrix (M z with z — 1,2, ...,n ), will 
depend on a certain number (j m ax(z)) of physical parameters each of which will have a particular distribution 
of values, i.e. an uncertainty, or fixed offset. These uncertainties on each parameter's value will have an effect 
on the Mueller matrix of the element that can be expressed by the product of the error on the parameter (5pj tZ ) 
and a matrix that tells us how this parameter affects each element of the main matrix, referred to as "weight" 
matrices (nxj )Z ) in the following. In this way we can express the Mueller matrix of any clement as the sum of a 
"main" matrix (M) and j mal (z) "error" matrices (Spj ■ mj): 



M z fa + 6pi, . . . ,pj + 5pj n0iX ( z )) = M(pi, • • • iPj) + Spx ■ mi H h 6pj ■ mj H h Sp jmax{z) ■ m jmax(z) . (5) 

Now the error matrices, (5pj, z ■ mj z ), can be described by the first order of the Taylor approximation of M 
with respect to each parameter pj tZ as showed by Keller & Snik (2009) and Tyo (2002) EED This is only valid 
if (1) the errors are small, compared to the value of the corresponding main matrix value, and (2) the errors 
are independent from each other. In some cases the error in certain parameters can be such that the first order 
of the Taylor expansion will not suffice to properly approximate the contribution to the main matrix and that 
higher orders have to be considered, but as a first step we will use the simplified expression. 

Now, to find the total Mueller matrix of a particular system composed of n elements we would have to, 
following the Stokes formalism, multiply the corresponding Mueller matrices: 

l 

Mtot = ]\ [ Mz (Pl.*> • ' ■ >Pj m ax(z),z) + tPl,* ' m l,» 4 1" tiPj ma *(z),Z • m j m ax(z),z] • (6) 

Z —n 

Taking into account the assumptions (1) and (2) made above, the resulting expression can be simplified as 
follows: 

1 1 3'ma.W /z + 1 \ / h=n \ 

M tot « n Mz + • jM k uj n Mk ■ ( ? ) 

z=n z=n j = l \k=l ) \k= Z -\ ) 

2.2 Implementation of Polarimetric Modulation 

The M&m's code is designed to simulate the polarimetric measurement process from calculating the Mueller 
matrix of a system to giving the final accuracy with which a certain Stokes parameter is measured. The procedure 
used to simulate the polarimetric measurement is schematically shown in Figure [l] and detailed below. 

Once we have the expression describing how the errors 5pj. z affect the main Mueller matrix of the system, 
we can calculate the total Mueller matrix for each modulation state (i). This enables us to build the modulation 
(O) and demodulation (D) matrices, including the error components, and eventually propagate the error to 
the final output of the polarimetric measurement process. The type of modulations implemented on the code 
are only spatial modulation, temporal modulation or a combination of both. Spectral modulation or channelled 
spectropolarimetry has a demodulation procedure based on Fourier transforms and therefore can't be represented 
by a linear operation. 

Continuing with the nomenclature presented in Section |2.1[ let us consider an optical system composed 
of n elements, each of which depends on a certain, and generally different, number of parameters [jPj, z , with 
j = 1, 2, ... , jmax(z) and z = 1, 2, . . . , n) and a modulation scheme composed of m modulation states. Note that 
we will then have an h — Y) r —j jmax{z) total number of parameters for the system, and j — 1,2, .. . , j max (z) 
inside each z optical element. Table [T] shows, for clarity, the definition and ranges of the indexes we will use in 
the following. 



Input 

Optical and polarimetric system layout 
>Modulation scheme: i = 1, 2, n 



> Total Mueller matrix for each modulation state i : 
Mi(Pi, P2,-) + 5p 1 -m il +8p 2 -m i2 +... 



^Modulation Matrix 

° = O maln +Sp 1 -o pl +5p 2 -o p2 +.. 



> Demodulation Matrix 

V"main "mainV "main 



^Response Matrix X+AX = D-0 
7 



Output 

S* in = (X+AX) Si„ 



Figure 1. Schematic view of the M&m's computational procedure. The program starts with the input given in the master 
script. Then, making use of the M&m's library, it computes the Mueller matrices for the n elements and the final Mueller 
matrix of the system, for each modulation state i. Next step is to build the modulation matrix taking the first row of each 
main and weight matrices to propagate the error on each parameter into the modulation process. Then the modulation 
matrix O is split into a main matrix O ma i n and the corresponding matrices Oj z for each parameter of the system. The 
demodulation matrix is either obtained computing the pseudo-inverse of O ma i n or specified by the user. The response 
matrix of the polarimeter X is then obtained as the product of the modulation and demodulation matrices and gives the 
relationship bet ween the measured and the incoming Stokes vectors. The matrix AX can then be used to perform the 
error budgeting 1^1121 

Table 1. Definition of indexes 



Optical clement 


z 


z = 1, 2, . . . , n 


Parameter in z 


3 


j =0,1,.. .Jmax(z) 


Modulation state 


i 


i = 1,2, ... , to 



For each modulation state i we would have a total Mueller matrix for the system: 



i=l, Mtot.i = Mi + Sp lt x ■ mi, 



5pj,z ■ m; z,l 



+ % mM {n),n ' m j ma5c (n),n,l 



i = 2, M to t,2 = M 2 + Spi A ■ 1111,1,2 H h 5p jtZ ■ m i z 2 H h 8pj matc (n),n • ™j maJC ( n ) ,n,2 



TO, 



M, 



Spi 1 ■ mil- 



8p 3 



S Pj 



c(n),n,: 



Each Mtot.i represents the behavior of the instrument for each modulation state. There will be only a few 
Mueller matrices that will change with (i), since is only the modulator (spatial or temporal) that changes state. 
The code will be able to distinguish between those ones and the "static" ones so the calculation of the total 
Mueller matrix can be shortened. 

Now if we have an incoming Stokes vector (Si n ) going through the system, for each modulation state we will 
get an intensity (ij) measured at the detector. These intensities form a vector at the end of the modulation 



cycle, I E 



(Ii,l2, ■ ■ ■ ,I m )- A linear combination of these to intensities will allow us to infer the state of 



polarization of the incoming light. This linear relation between the m intensities measured by the detector and 
the incoming stokes vector is given by the so-called modulation matrix (O), 



Imeas — OS; n . (8) 

The modulation matrix O is an m x 4 matrix built by compiling the first row of every Mtot,i matrix. This 
can be understood remembering that it is in this row of the Mueller matrix where we find how the components of 
the incoming Stokes vector get transformed into output intensity Eq. [5] Given the structure that these Mueller 
matrices now have, i.e. being split in a sum of matrices, we will obtain a modulation matrix that is also a sum of 
matrices. We will get a main modulation matrix (O ma i n ) obtained by compiling the first rows of main matrices 
Mi and h matrices Oj z , obtained by compiling the first rows of weight matrices m.- }Z i 1 each one associated to 
one physical parameter of the system. 

O = O main + Spi iX ■ oj.,1 H h 5p jtZ ■ Oj z H h 5p jmax ( n ) !n ■ o jmax ( n ) in . (9) 

Omain tells us how the intensity is measured without taking into account the deviations of the parameters. 
The "separated" propagation of the errors to this point allows us to analyze how they affect the intensity 
measurement directly, so we can easily detect the elements and parameters that have more impact on the 
measurement before the demodulation process. 

It is at this point of the process when we can implement detector errors such as bias drift and nonlinearities 
(see Section [3] for further description). 

Once we have the measured intensity vector the demodulation process starts. The basic goal of this process 
is to infer S- ln from the set of intensity measurements I mea s- It is clear from Eq|8] that this can be done by 
the inverted matrix of O, which we will refer to as D, the demodulation matrix. This is a trivial mathematical 
problem only if O is a square matrix, which is the case for modulation schemes composed of four measurements. 
If the modulation scheme comprises a larger number of modulation states, i.e. m > 4, the solution for the 
inversion of O is no longer unique and we have an infinite number of demodulation matrices that fulfill the 
condition DO = 1. Then the new problem is to find the demodulation matrix that optimizes the efficiencies of 
the systemHEH In 2000, del Toro Iniesta & Collados^ showed that the solution for this optimization problem is 
found by means of the Moore-Penrose pseudo-inverse matrix: 

D = (O t O)- 1 O t . (10) 

At this point we need to make some considerations. Given expressions j9j and ^ we would now have to 
calculate the Moore-Penrose pseudo- inverse matrix of a sum of a large number of matrices, depending on the 
total number of parameters in our system. Furthermore, we want to keep the errors 5pj yZ as unknown variables 
up to the end of the process so we can analyze the system in a generic way, without considering specific values 
for the errors. There are various mathematical formulae and theorems to obtain the inverse of a sum of matrices 
but , as far as we know, there is no way to obtain analytically the inverse of a structure like this particular 
one, specially with a variable number of matrices. The problem of how errors propagate through the non-linear 
process of matrix inversion has been analyzed by Asensio Ramos & Collados in 2008 for gaussian errors^ but 
systematic errors have a complete different distribution which makes them unsuitable for this procedure. A 
simplified approach is to obtain D by means of an "ideal" or "best known" modulation matrix, which would be 

Omain m OUT Case. 

This can be justified arguing that during data reduction one has to choose a demodulation matrix. The best 
knowledge one has of the modulation matrix is obtained after modeling (or calibrating). This matrix would be 
O ma i n for us, since it represents de "ideal" behavior of the modulation matrix. In this way we can also propagate 
the error in the parameters as unknown variables through the demodulation process. 

Therefore: 

D* = fO T ■ O ■ ) _1 O t • (11) 



There is also the possibility of using a pre-defined demodulation matrix which, even though it may not be 
the optimal one, will be required in many practical situations. This can be easily implemented in the code by 
providing the demodulation matrix together with the input and not deriving it from the calculated modulation 
matrix. Also the code needs to be suitable for partial polarimeters, i.e. polarimeters only measuring a part of the 
Stokes vector, for which the modulation matrix when written as a m x 4 matrix is singular and non-invertible. 
Therefore the code would have to be able to reduce the dimensionality. 

Once we have the demodulation matrix we can finally find the vector that will represent the "measured" 
incoming Stokes vector, Sf n , 

s; n = DI meas . (12) 

Note that this vector is not a Stokes vector since is obtained through the matrices O and D which are not 
Mueller matrices. 

Ideally this output will be exactly the Stokes vector that comes in the polarimeter so that: 

= iS in . (13) 

Of course this is not the case since no instrument is exempt from creating some uncalibrated instrumental 
polarization and cross-talk along the measurement process. To take this into account it is useful to introduce 
the concept of a response matrix This matrix, which is not a Mueller matrix, has dimensions of 4 x 4 and 
can be defined as the one representing the measurement process i.e. it relates the incoming Stokes vector (Si n ) 
with the vector S; n that we obtain in the end as the "measured" Stokes vector 

S? n = XS in . (14) 

Again, in ideal conditions X would be the identity matrix, but since we have errors in the process it will 
have values different from unity. In the end, what we aim to do with the M&m's code is to model this response 
matrix, so find the range of values its elements can take on, and this is done by propagating our error-dependent 



matrices Oj jn up to the end of the inversion process. Considering equations |8j [12] and 14 it is obvious that we 
can obtain X by means of D and O as follows: 

X = DO. (15) 



If we now introduce the expression for O, Eq. 15 becomes 



X = D[O main + <5p M • o 14 H h 5p j>z ■ o j:Z H h <% max („),„ ■ Oj m<uc ( n ), n ] . (16) 

which applying the distributive property gives: 

X = DO main + 5p ia ■ (Do M ) H h Sp j>z ■ (Doj, z ) H h <% max („),„ ■ (Do jmax(n ), n ) . (17) 

If the demodulation matrix is the pseudo-inverse of the main modulation matrix, i.e. we have computed D 
instead of using a pre-defined one, the first term of the right side of the equation is the identity matrix and each 
(Doj z ) term represents a variation of the response matrix depending on the corresponding Pj, 2 . 

X = 1 + Sp hl ■ AXi,! + • ■ ■ + 6p jtZ • AXj, z + • • • + <%_(„),„ • AX jmax(n) ,„ = 1 + AX (18) 



In this way, we have propagated the errors up to the final matrix that models the complete measurement 
process. 



2.3 Error budgeting 



The matrices AXj z and AX, from Equation 18 can be seen as "accuracy matrices" to be directly used for 
error budgeting by comparing them with the requirements set for the polarimeter. 

This can be done in different ways such as comparing the complete accuracy matrix (AX) with a "require- 
ments matrix" that one builds from the scientific requirements^ (Eq. 19): 



AX < AX r 



(19) 



can be considered. 



To simplify the process only the most stringent components of the AX req 

Another approach is to assume an input Stokes vector (Si n! def 
between the this input and the measured Stokes vectors (S? n ) obtained by multiplying with AX (Eq 



and compute the difference or vector norm 

"TBI. 



20 



4 

£(AX.S in , def )^<lkH 2 - (20) 

k=l 

The requirement on polarimetric accuracy is then described by the vector norm ||e||, and should be fulfilled 
for all possible input Stokes vector. 

2.4 Implementation of Calibration 

The code can be set up such that it describes the full procedure of calibration, after the Mueller matrices for the 
calibration optics including error terms have been introduced. These results should then be fed back into the 
code to describe the calibrated polarimetric measurements. However, the error propagation from the calibration 
to the measurement is not straightforward and has only be mathematically described in the case of Gaussian 
errors.^ Possibly a complete Monte Carlo simulation has to be performed to fully model this process. 

In any case, real calibration results and the estimated errors thereupon can be implemented directly by 
replacing the pertinent part of the modeled Mueller matrix train. 



3. OVERVIEW OF POLARIMETRIC ERRORS: THE PHYSICS 

Now that we have defined the mathematical procedure we need to introduce the physics behind the error 
modeling, i.e. determine mj z . An error budgeting code like M&m's is only useful once it is able to consider 
almost any kind of error to any kind of optical element. Based on the current literature, an extensive library is 
therefore created out of which most relevant Mueller matrices and error matrices can be taken or computed. 

Errors affecting polarimetric measurements can be widely ranging in both type and scale and can arise from 
many different sources like missalignements, manufacturing errors or variations of the working temperature. 
In polarimetry, every element before the polarimeter, even a simple glass plate, can affect the measurement. 
Therefore a careful analysis of each element must be done in order to model its behavior as accurately as 
possible. 

There are two main ways in which this can be done. One way is based on the direct characterization of 
the elements in the lab, or what could be seen as the empirical approach. This method, of course, provides a 
quite complete knowledge of the particular element one is considering but lacks of potential for generalization 
and it is limited by the accuracy of the measuring system. Another approach is to start from the general 
analytical expression of the element's Mueller matrix depending on various physical parameters and derive the 
error matrices as a first order Taylor expansion to each parameter. In our case we consider both possibilities, so 
the code can estimate an approximation of the weight matrices for each parameter and it can also include error 
matrices that are already determined. 



3.1 Parameters 

Parameters can be separated in two groups, the ones affecting the Mueller matrix of only one particular clement, 
e.g. the birefringence of the material of a retarder or the extinction ratio of a polarizer, and those which are 
common for all the system, e.g. wavelength and temperature. The latter will be represented in the code by global 
parameters. This means that all Mueller matrices inasmuch they depend on wavelength and temperature will 
automatically have error matrices associated to them that scale with wavelength or temperature. Thus the varia- 
tion of system performance with wavelength and temperature can be estimated from the code's output. It is clear 
that with only a first-order approximation of Mueller matrix variation with wavelength, only relatively narrow 
band-widths can be considered. If one element has a different temperature than the others, its Mueller matrix 
at that temperature and corresponding error matrix for temperature changes can be also inserted separately. 

3.2 Library of Mueller Matrices and Error Matrices 

The following describes the implementation of Mueller matrices of common optical elements. It also describes 
which error sources can be used in the code and under which assumptions. 

3.2.1 Rotations 

This is a special case because it is not an element per se. A single rotation Mueller matrix describes the rotation of 
the [Q, U] coordinate system. The only error in this process is an offset to this rotation (or a certain distribution 
of offsets). For the case of a freely rotating element, two rotation matrices need to be inserted, one before and 
one after the element (which can be a compound of several Mueller and error matrices). The amount of rotation 
and the error thereof is identical, but have an opposite handedness, unless the elements in between modify Stokes 
coordinate system like an odd number of mirrors do. 

Rotation is often used in temporal modulation, so the variation of the Mueller matrix with modulation state 
is generically implemented in this function. 

3.2.2 Linear Polarizers 

The ideal Mueller matrix of a linear polarizer assumes that linear polarization in the +Q direction gets trans- 
mitted. In the case of spatial modulation by using a polarizing beam-splitter, this direction can be switched to 
— Q in the input script that handles the modulation sequence. 

The following error terms are implemented: 

• Limited extinction ratio. The extinction ratio is here defined such that its value is <1. The error matrix 
for limited extinction ratio is described by Snik (2006) This error necessarily has a first and second 
order. 

• Depolarization cf. Nee et al. (1998)PI 

• Field-of-view effects. This is described with the 3D geometry of the polarizers axes compared to the input 
polarization. 

3.2.3 Retarders 

The most frequently used retarders are based on birefringent (liquid) crystals. The retardance can be computed 
once the material's birefringence and the plate's thickness is known: 

5 = 2n {ne -^ ) - d . (21) 
A 

The refractive indices for many crystals are listed by Ghosh (1998)P 
Various error sources are readily derived from this: 



• Variation with wavelength. 



Variation of the thickness. 



• Offset of the birefringence. 

• The variation of the birefringence with temperature can be computed from thermo-optic coefficients listed 
by Ghosh (1998). See Hale & Day (1988). 8 Also the thermal expansion of the plate is taken into account. 

• The variation of retardance upon a change of incidence angle and azimuth, as described by the equations 
in Evans (1949)P 

• Dichroism, as determined by the Fresnel equations for the two polarization directions along the two crystal 
axes. 

• Polarized fringes due to multiple reflections within the bircfringent plate manifest themselves as a sinusoidal 
partial polarization pattern as a function of wavelength, as well as a wavelength-dependent modification 
of the retardance. These effects are calculated by using the Berreman calculus as described by Weenink et 
al. (2011)P^ The mitigation of these effects by using AR coatings can be implemented ad-hoc. 

Since most wave-plates used for instrumentation are actually compounds of two or more layers of different 
materials, it makes sense then to consider compound wave plates as unit elements in the system. Several 
frequently used wave plates will be predefined in the M&m's library: achromatic and superachromatic wave plates 
and liquid crystals like Liquid Crystal Variable Retarders (LCVRs) en Ferroelectric Liquid Crystals (FLCs). For 
the first cases, the alignment between the two plates becomes an important error parameter. For the liquid 
crystals, the switching angles (that often drive the modulation sequence) will be subject to errors. 

Also a standard Fresnel rhomb will be implemented. The main error sources there will be: 

• Variation of retardance with incidence angle and azimuth. 

• Stress birefringence. 

3.2.4 Mirrors 

The Mueller matrix of a mirror is based on the models applied to accurate ellipsometry performed by Van Harten 
et al. (2009) According to it, a mirror can be characterized knowing the following parameters: the incidence 
angle (a) and wavelength (A) of the incoming light, the complex index of refraction of the metal, and, if present, 
the thickness and index of refraction of the dielectric film layer on top of the mirror's surface. This layer can be 
an artificial overcoating or can occur naturally due to the growth of the aluminum oxide layerP^O The equations 
that lead to a mirrors Mueller matrix are numerous and complex and directly describe the dependence of the 
Mueller matrix on various physical parameters: 

• Variation with wavelength due to the variation of the refractive indices. 

• Variation with incidence angle and azimuth. 

• Aging which involves a growth of the layer of dielectric material on the mirror. Also a dust layer can be 
described with an effective growth of this layer, see Snik et al. (2011)p31 

3.2.5 Glass Components 

It is often assumed that glass components like lenses do not modify or create polarization, and therefore the 
ideal Mueller matrix is the unit matrix. Various error sources can however exist: 

• Glass components can introduce linear polarization if a ray of light hits it at a non-normal incidence. This 
is fully determined by the Fresnel equations, after linearizing those for the incidence angle. It needs to be 
computed at every glass-air interface. 

• Stress birefringence occurs in every piece of glass. A certain direction (or distribution) of the stress tensor 
needs to be assumed. 



3.2.6 Detectors 

Although the action of a detector cannot be described by a Mueller matrix, its non-ideal properties can adversely 
affect the polarimctric accuracy. Fortunately, these error terms can easily be implemented using a similar 
formalism. 

It is assumed that an ideal detector linearly converts the intensity to a data number. This corresponds with 
a fully multiplicative term comparable with a diagonal Mueller matrix. Two important error terms are: 

• Bias drift or uncorrected stray light. 

• Detector non-linearity This obviously needs to be described with a second order (or higher) error term. 
Keller (1996)P^ has shown that such errors can couple with instrumental polarization to create spurious signals. 



4. OVERVIEW OF THE CODE 

The M&m's is a code written in Python which makes it accessible to a broad audience cost-free. The name arises 
from the basic structure of the Mueller matrices (Eq. [5]) that is propagated all along: the "main" matrix (M) and 
the "weight" matrices (m).It is structured in two main parts, the "input/master script" and the u M&m's library". 



The input/master script represents the input to the program. In it, the user has to specify all the information 
about the set up to be simulated, e.g. type and order of the optical elements on the light path, modulation states, 
global parameters, rotations and misalignments between them, etc.... In this input one can also provide particular 
matrices to be used instead of those from the library such as specific calibrated error matrices for a certain element 
or demodulation matrices. The user also needs to specify the mode in which the code should operate. 

The M&m's library is a compilation of functions that can be separated in two groups according to their 
functionality. 

The first group is composed of "element functions" computing the Mueller matrix of typical optical elements 
in terms of the physical parameters in which the element depends on and deriving the corresponding weight 
matrices. The library contains functions to model mirrors, glass elements, waveplates and polarizers among 
others. 

The second group is formed by the "mode functions" which, depending on the mode selected, operate with 
the Mueller matrices obtained for the elements in order to get the desired output. Figure [2] shows a flow diagram 
of the information through the different elements and stages of the code. 



Input/master script 



> Optical and polarimetric 
system layout 

> Modulation scheme: 
i = 1, 2, m 

> Operational modes: 

• ModelMuellerMatrix 

• Mode/Measurement 

• Mode/Calibration 



M&m's library 



> Element funct. 



Mode funct. 



Output 



> ModelMuellerMatrix 
M,„,(Pi, p 2 , - ) + <5p 1 -m 1 +<5p 2 -m 2 +... 



ModelStokesMeasurement 
S* rn = (X+AX) Si n 



> ModelSystemCalibration 
O cal = (O main )«'+ 6p 1 {o pl ) al +... 



Figure 2. Flow diagram of the M&m's operation. The master script sends the input information to the library's selected 
element functions to be used. The output of this operation is sent to the master script which then sends the matrices to 
the selected mode function producing the final required output. 



With the aim of making it suitable for different types of analysis it operates in three modes: 



• ModelMuellerMatrix 

• ModelStokesMeasurement 

• ModelSystemCalibration 

The first mode, ModelMuellerMatrix, computes the total Mueller matrix of a given system. It takes the 
matrix for each element that come as the output of the element functions and then multiplies the matrices 
following the formalism explained above. It gives the total Mueller matrix for the system as an output with the 
structure shown in Eq. [7j This is useful to analyze how the different element/parameters will affect the final 
behavior of the polarimetric set up but also to estimate the Mueller matrix of optical systems upstream, i.e. the 
telescope in front of the polarimeter. 

The second mode, ModelStokesMeasurement, computes the vector that it is measured at the end of the process 
(S* ut ), i.e. it propagates the error matrices through the modulation and demodulation operations. In the end we 
get the correspondence between the incoming and the measured Stokes vector by means of the response matrix 
X which contains the effect of the physical parameters of the system on the measurement, represented by the 
propagated weight matrices (AX). 

The third mode, ModelSystemCalibration computes the polarimetric modulation matrix as obtained after 
calibration and the errors pertaining to these results. This mode has not yet been implemented since, as we 
pointed out before, the error propagation through this operation to a calibrated measurement model needs to 
be investigated. 

5. DISCUSSION AND OUTLOOK 

In this manuscript we have presented the mathematical, physical and operational basis of the M&m's code for 
simulating the performance of a given polarimetric system. The mathematical framework developed by Keller & 
Snik allows us to propagate the contribution of the different errors separately up to the end of the measurement 
process. The main advantage of this approach is that it allows us to, at any point of the measurement process, 
identify the critical sources of deviation in the measurement without having to specify any particular error for 
the parameters. Furthermore the matrix (AX), which is dependent on the Spj scalar errors, can be directly 
compared to the accuracy matrix required for the system to set the tolerances on the parameters, because we 
have propagated Spj as unknown variables through the process. 

The main disadvantage is that it relies on assumptions (1) and (2) presented in Sec. |2.1| (errors are small and 
independent) when aproximating the weight matrices by the first term of the Taylor expansion of the Mueller 
matrix. Higher order terms may therefore need to be implemented. 

The code will be verified by comparing its results with those from a full Monte Carlo simulation (which can 
be done by using the same library of Mueller matrices), and with the results of a dedicated lab set up. At the end 
only a Monte Carlo simulation or lab measurements can fully characterize the system, but the M&m's approach 
provides a better insight on the system's dependencies. This allows us to easily detect, early in the design 
process, which ones will be the limiting elements/parameters when analyzing the accuracy that different designs, 
considered for a polarimetric system, can provide. It also can easily estimate the instrumental polarization 
introduced by optical systems, keeping also track of which elements have more impact. 

The code aims to be open, so it can be used by anyone, general, so it can be applied to any polarimetric 
design, and complete in terms of error modelling. In this sense, we will keep the library open so it can be easily 
updated with new matrices for elements and errors. 

The code is still in an early stage of development and it will be made accessible to the public as soon as it 
has been verified, and most relevant errors have been implemented. 

The following upgrades to the code are still under consideration: 

• A complete description of a distribution (random or systematic) of the physical parameters that drive the 
errors. This way, after propagation, the error contributions can be added in an RSS fashion, see Keller & 
Snik (2009). 



• Adding the possibility of introducing higher-order error terms. 

• Full implementation of the simulation of the calibration process and the calibrated measurement process. 
It may be necessary to implement a Monte Carlo simulation of the calibration process to fully propagate 
the systematic errors, as currently only Gaussian errors can be propag atecP . 

• Integration with the optical design, such that all Mueller matrices are explicitly dependent on the local 
incidence angles and azimuths. For converging beams, an average Mueller matrix plus a distribution of 
values around that can thus be obtained. A link with ZEMAX would be a logical choice for this. A simple 
simple solution would be to parameterize the distance to the pupil plane for each element and supply a 
range of incidence angle and azimuths. 
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